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Abstract 

In this paper we present the theoretical framework needed to justify the use of 
a kernel-based collocation method (mcshfrec approximation method) to estimate the 
solution of high-dimensional stochastic partial differential equations (SPDEs). Using 
an implicit time stepping scheme, we transform stochastic parabolic equations into 
stochastic elliptic equations. Our main attention is concentrated on the numerical 
solution of the elliptic equations at each time step. The estimator of the solution of the 
elliptic equations is given as a linear combination of reproducing kernels derived from 
the differential and boundary operators of the SPDE centered at collocation points to 
be chosen by the user. The random expansion coefficients are computed by solving a 
random system of linear equations. Numerical experiments demonstrate the feasibility 
of the method. 
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1 Introduction 



Stochastic partial differential equations (SPDEs) frequently arise from applications in areas 
such as physics, engineering and finance. However, in many cases it is difficult to derive an 
explicit form of their solution. Moreover, current numerical algorithms often show limited 
success for high-dimensional problems or in complex domains - even for deterministic partial 
differential equations. The kernel-based approximation method (meshfree approximation 
method [4, 11, 21]) is a relatively new numerical tool for the solutions of high-dimensional 
problems. In this paper we apply - to our knowledge for the first time - such a kernel-based 
collocation method to construct numerical estimators for stochastic partial differential equa- 
tions. Galerkin-type approximation methods based on the eigenvalues and eigenfunctions 
of the underlying differential operator are currently very popular for the numerical solution 
of SPDEs [7, 17, 14]. With the kernel-based meshfree collocation method introduced here 
explicit knowledge of these eigenvalues and eigenfunctions is not required since the kernels 
can be directly obtained as Green kernels of the differential operators [10, 9]. Stochastic 
collocation methods using a polynomial basis are also frequently found in the literature 
[2, 18]. These methods usually require the collocation points to lie on a regular grid. In our 
method the collocation points can be placed at rather arbitrarily scattered locations. This 
allows for the use of either deterministic or random designs such as, e.g., uniform or SoboP 
points, but also for adaptively chosen locations. In this paper we do not study the design 
aspect of our algorithm, but reserve this important aspect for future work. Another advan- 
tage of using a meshfree method is its ability - also not explicitly pursued here - to deal 
with problems on a complex domain T> C M. d , d > 1, by using appropriately placed colloca- 
tion points. Another advantage of this method is its high efficiency, in the sense that once 
certain matrices are inverted and factored we can compute, essentially for free, the value of 
the approximated solution at any point in the spatial domain and at any event from sample 
space. In particular the method is suitable for simulation of a large number of sample paths 
of the solution. In this article we present only a general framework for this new numerical 
method and prove weak convergence of the corresponding schemes. We conclude the paper 
with a numerical implementation of this method applied to a one-dimensional stochastic 
heat equation with Dirichlet boundary conditions driven by an additive space-time white 
noise (colored in space). Much more details, as well as some of the aspects just mentioned, 
will be discussed in Qi Ye's Ph.D. thesis [23] or in future publications. 

1.1 The method in a nutshell 

Assume that T> is a regular open bounded domain in M. d (see Appendix B), and let Ti 
be a separable Hilbert space of functions defined on V. Also, let (£lw, Fwi i^t}, Pw) be 
a stochastic basis with the usual assumptions. We consider the following parabolic Ito 
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equation 

' dU t = AU t dt + adW t , in V, < t < T, 
BU t = 0, on &D, 

U = u Q , 

where A is a linear elliptic operator in T~L , B is a boundary operator for Dirichlet or Neumann 
boundary conditions, uq G H, and W is a Wiener process in H, with mean zero and spatial 
covariance function R given by K(W(t,x)W(s,y)) = min{t, s}R(x, y), x,y G V, t, s > 0, 
and a > (see for instance [5]). 

We assume that equation (1.1) has a unique solution U G L2(^ty x (0, T);H). 

The proposed numerical method for solving a general SPDE of the form (1.1) can be 
described as follows: 

51) Discretize (1.1) in time by the implicit Euler scheme at equally spaced time points 
= t o <h< ... <t n = T, 

U tj - U ti _ x = AU tj 5t + aSWj, j = 1, . . . , n, (1.2) 

where St := tj — and SWj := W tj — W^^. 

52) Since it follows from (1.2) and the definition of Brownian motion that the noise in- 
crement SWj at each time instance tj is independent from the solution Ut j _ 1 at the 
previous step, we simulate the Gaussian field with covariance structure R(x, y) at a 
finite collection of predetermined collocation points 

X v := {xi,--- ,x N } C V, X dT> := {x N+ i,--- ,x N+M } C &D. 

53) Let the differential operator P := I — StA, and the noise term £ := aSWj. Since £ is a 
Gaussian field with E(£a.) = and Cov^a-^y) = o~ 2 StR(x,y), equation (1.2) together 
with the corresponding boundary condition becomes an elliptic SPDE of the form 

/ + £, inP, 

(1.3) 

0, on &D, 

where u := is seen as an unknown part and / := Ui-_ x and £ are viewed as given 
parts. We will solve for u using a kernel-based collocation method written as 

N M 

* * 

u(x)ttu(x):=^2c k P 2 K(x,x k ) + ^2c N+k B 2 K(x,x N+k ), x G V, (1.4) 

k=l k=l 

* * * 

where K is a reproducing kernel and the integral-type kernels K, P2K, B2K are defined 
in Lemmas 2.1 and B.2. The unknown random coefficients c := (ci, • • • ,c/v+m) T are 
obtained by solving a random system of linear equations (with constant deterministic 
system matrix and different random right-hand side) at each time step. Details are 
provided in Section 3. 
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S4) Repeat S2 and S3 for all j = 1, . . . , n. 



Obviously, many other - potentially better - time stepping schemes could be applied 
here. However, as mentioned earlier, we focus mainly on step S3 and are for the time 
being content with using the implicit Euler scheme. Naturally, the rate of convergence of 
the above numerical scheme depends on the size of the time step St, and the fill distance 
hx := supa.gx' minajj.gXpuXax, \\x — Xk\\2 of the collocation points. We support this statement 
empirically in Section 4. We should mention that even for deterministic time-dependent 
PDEs to find the exact rates of convergence of kernel-based methods is a delicate and 
nontrivial question, only recently solved in [12]. We will address this question in the case 
of SPDEs in future works. 

The fundamental building block of our mesh-free method is the reproducing kernel 
K : D x D — > R and its reproducing- kernel Hilbert space Hic(P) (see Appendix A for 
more details). By the very nature of such a kernel-based method, the approximate solution 
lit- , j = 1, . . . ,n, must live in Hr-('D). Thus, we make the following standing assumption 
throughout the paper: 

The solution U of the original equation (1.1) lies in a Hilbert space T~L which can be 
embedded in the reproducing kernel Hilbert space H/<(D). 

Usually, % is a Sobolev space % m (D), for some positive m. In this case it is possible to 
choose an appropriate kernel K such that the above embedding assumption is satisfied. For 
a general discussion of existence and uniqueness of the solution of problem (1.1) see, e.g., 
[19, 6, 5]. 

2 Reproducing-kernel collocation method for Gaussian pro- 
cesses 

In this section we briefly review the standard kernel-based approximation method for high- 
dimensional interpolation problems. However, since we will later be interested in solving a 
stochastic PDE, we present the following material mostly from the stochastic point of view. 
For further discussion of this method we refer the reader to the recent survey papers [20, 8] 
and references therein. 

Assume that the function space Hk('D) is a reproducing-kernel Hilbert space with 
norm || • \\k,V and its reproducing kernel K £ C(T> x V) is symmetric positive defi- 
nite (see Appendix A.l). Given the data values {yj}j = i C R at the collocation points 
Xt> '■= {xj}^ =l C P of an unknown function u £ H^-(P), i.e., 

y j = u(x j ), xj = (xij,--- ,x dJ ) e V C R d , j = l,...,N, 

the goal is to find an optimal estimator from H^-(P) that interpolates these data. 
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Definition 2.1 ([3, Definition 3.28]). A stochastic process S : T> x Vt — > M is said to be 
Gaussian with mean fx : D — > R and covariance kernel f:DxP-*lona probability 
space (f^J 7 , P) if, for any pairwise distinct points Xp := {x\,--- ,xn} C T), the random 
vector S := (5' £Cl , • • • , S XN ) T is a multi-normal random variable on (Q, J 7 , P) with mean /x 
and covariance matrix 0, i.e., S ~ Af(/J,,<b), where /x := (//(a;i),--- ,[i(xn)) t and <t> := 
(*(ajj,a!fc))^.= r 

2.1 Data fitting problems via deterministic interpolation and simple krig- 
ing 

In the deterministic formulation of kernel interpolation we solve an optimization problem 
by minimizing the reproducing-kernel norm subject to interpolation constraints, i.e., 

u K = argmin {||u||k,x> s.t. u(xj) = yj, j = 1, . . . , N} . 

In this case, the minimum norm interpolant (also called the collocation solution) uk{x) is 
a linear combination of "shifts" of the reproducing kernel K, 

N 

u K {x) :=^2c k K(x,x k ), xeV, (2.1) 

k=l 

where the coefficients c := (ci,--- ,cn) t are obtained by solving the following system of 
linear equations 

Kc = y , (2.2) 

with K := (Kix^Xk))^^ and y := (yi,--- ,?/at) t . 

For simple kriging, i.e., in the stochastic formulation, we let S be a Gaussian process 
with mean and covariance kernel K on some probability space (f2, J-, P). Kriging is based 
on the modeling assumption that u is a realization of the Gaussian field S. The data values 
2/1, ... , are then realizations of the random variables S X1 , . . . , S XN . The optimal unbiased 
predictor of S x based on S is equal to 

N 

U x := ^ c k {x)S Xk = argmin E\U - S x \ 2 , 

fc=l C/espan{5 a!j }^ =1 

where the coefficients c(x) := (ci(x), ■ ■ ■ ,cn(x)) t are given by 

c{x) = K- l k(x) 

with k(x) := (K(x, x%), ■ ■ ■ ,K(x,xn)) T and the same matrix K as above. We can also 
compute that 

^(U x \S Xl = yi, ■ ■ ■ , S XN = y N ) = u K {x). 

Note that, in the kriging approach we consider only the values of the stochastic process S 
at the collocation points, and view the obtained vector as a random variable. However, if we 
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view S as a real function, then P(S* G Etft-('D)) = by [16, Theorem 7.3]. A simple example 
for this fact is given by the scalar Brownian motion defined in the domain T> := (0, 1) (see, 
e.g., [10, Example 5.1]). This means that it is difficult to apply the kriging formulation 
to PDE problems. Next we will introduce a new stochastic data fitting approach that will 
subsequently allow us to perform kernel-based collocation for stochastic PDEs. 

2.2 Data fitting problems via a new stochastic approach 

From now on we will view the separable reproducing- kernel Hilbert space (T>) as a sample 
space and its Borel a-field ^(Hft-(D)) as a cr-algebra to set up the probability spaces so that 
the stochastic process S x (u) := oo{x) is Gaussian. We use the techniques of [13, 16] to verify 
Lemma 2.1, which is a restatement of [16, Theorem 7.2]. This theoretical result is a general- 
ized form of Wiener measure defined on the measurable space (C[0, oo), ^(C[0, oo))), called 
canonical space, such that the coordinate mapping process W x (u) := oo{x) is a Brownian 
motion (see, for instance, [15], Chapter 2). 

Lemma 2.1. Let the positive definite kernel K G C(V x V) be the reproducing kernel of 
the reproducing-kernel Hilbert space ^(P). Given a function /i G Hk^) there exists 

a probability measure P M defined on {VLk-,^Fk) '■= (Hj^p), ^(H^(D))) such that S x (uj) := 

* 

u(x) is Gaussian with mean [i and covariance kernel K on (O^, J 7 ^:, P M ), where the integral- 
type kernel K of K is given by 

K(x,y):= / K(x,z)K(y,z)dz, x,y£T>. 
Jv 

Moreover, the process S has the following expansion 

oo 

S X = Y^ CkVXke k (x), x G V, F^-a.s., 
k=i 

where {A^}^^ an ^ { e fc}fcLi are ^ ne eigenvalues and eigenf unctions of the reproducing kernel 
K, and £fc are independent Gaussian random variables with mean fi^ := (/i, \/Afc e fe)-ftT,X' an d 
variance A^, k G N. 

Before we prove Lemma 2.1 we remark that we have introduced the integral-type kernel 
K for convenience only. As seen later, in order to "match the spaces", any other kernel 
that "dominates" K (in the sense of [16]) could play the role of the integral- type kernel K. 

Proof. We first consider the case when \x = 0. There exist countably many independent 
standard normal random variables {£fc}fcLi on a probability space (0^, J^, P^), i.e., ~ 
i.i.d. A/"(0, 1), k G N. Let {Afc}-^ and {ek}^ =1 be the eigenvalues and eigenfunctions 
of the reproducing kernel K as in Theorem A.l. We define S := X]*^=i£fc^fc e fc P^-a.s. 

Note that S is Gaussian with mean and covariance kernel K. Since E(Sfc^=i Cfe-^fc) — 
SfcLi Var(^)Afe = J2T=i < oo indicates that YlkLi |£fc\/^fc| 2 < 00 P^-a.s., Theorem A.l 
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shows that S(-,lo) E Hx(f) P^-a.s. Therefore S is a measurable map from (0^,7-"^) into 
(Qk,J~k) by [3, Chapter 4.3.1] and [16, Lemma 2.1]. On the other hand, the probability 
measure P° := P^ o 5" 1 (also called the image measure of P^ under S) is well defined on 
(0,K, Fk), i-e., P°(^4) := P^(S ,_1 (A)) for each A E J 7 ^. Hence, S is also a Gaussian process 

with mean and covariance kernel K on (Qk, J~k, P )- 

Let S 11 := S + ii on (^.F^P ). Then E(S£) = E(S X ) + /x(a;) and Cov(S£,S(;) = 
Cav(S x , S y ) with respect to P°. We define a new probability measure P M by P At (A) := 
P°( J 4 — /i) for each A E J 7 ^-. It is easy to check that Q,k + [i = Hft-(f]) = Q,k and 
{/i + A : A E J~k\ = ^(H/f (fi)) = J 7 ^-. Thus S" is Gaussian with mean \i and covariance 

kernel K on (0^-, J 7 ^, P^). 

Moreover, since \x E H^(fi), it can be expanded in the form [i = YlkLi Afc\/^A; e fc> where 
Afc := (m, yfhe^Kp, so that = YlkLiif^ + V / AfeCfc)\Afce fe . But then Cfc ~ Afe + \Afc6 ~ 
■^(/U/o-^fc) are independent on (f^, .Fx, P M ). □ 

According to [3, Theorem 4.91], we can also verify that the random variable Vf(co) := 
(w, /)k,v, f £ Hr-('D), is a scalar normal variable on (Sl^, Fr-, P M ), i.e., 

~jV(m/,aJ), weO x = H K (D), 

where mj := {^,J)k,v and cry := ||/||l 2 (©)- Therefore the probability measure P M defined 
in Lemma 2.1 is Gaussian. 

Let : M N — > R be the joint probability density function of S^, • • • , defined on 
(f^K, Fr, P^ 1 ). Then it is a normal density function with mean /j, := (fj,(xi), ■ ■ ■ ,fi(xj\f)) T 

and covariance matrix K := (K(xj, scj;))^ =1 . In analogy to the kriging formulation we can 
find the optimal mean function fi E H^(P) fitting the data values y := (yx, • • • , vn) T , i-e., 

ft:=*k T K- l y = sup p£(y ) = sup F»(S = y ), 

* * * 

where := (i^(a;, x\), • • • , X(o;, xn)) T - 

We now fix any x £ T>. Straightforward calculation shows that the random variable S x , 

given S Xl , • • • , S XN , defined on (f^, Fr, P^) has a conditional probability density function 

1 {v-m£(v))»\ _ n _ mJV 



: = / s — exp I , vEM, t)£ 

where m x (v) := //(a;) + fe(a;)- r K~ 1 (?; — m' 1 ), m/ 4 := (/i(a;i),--- ,ii{x^)) T , and <r(a:) 2 
K(x,x) — k(x) T K^ 1 k(x) . Then the optimal estimator that maximizes the probability 

maxP M ({w E 0,r '■ oj(x) = v s.t. oj(x{) = y±, ■ ■ ■ ,uj(xn) = Vn}) 
= maxP /i (S x = v\S Xl =yi,--' ,S XN = Vn) 
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is given by 

u(x) := y{x) = argmax p x (v\y Q ). 
Proposition 2.2. With the above notations, the following equality holds true 

J>£(«(aOl!/o) = SU P P£Hlto)- (2.3) 

Moreover, for any e > 0, 

sup P^(|fi(a;)-u(aO|>e)< sup ¥^£ x ) = erfc ( \ ) , (2.4) 
ii&t K {V) A*eH x (x>) \\J2o-(x)J 

where 

£ e x := {u £tl K ■■ \u{x) - u(x)\ > e s.t. u(xi) = y 1 , ■ ■ ■ ,u>(x N ) = y N } . 

Identity (2.3) follows by direct evaluations. Consequently, taking into account that S is 
Gaussian, inequality (2.4) follows also immediately. 

Remark 2.1. Instead of giving a deterministic (or strong) error bound for the proposed 
numerical scheme, we provide a weak type convergence of the approximated solution u to 
the true solution u, as stated in Proposition 2.2. In fact, inequality (2.4) can be seen as a 
confidence interval for the estimator u with respect to the probability measure P^ 1 . 

In the next section we generalize this stochastic approach to solve elliptic PDEs. 

3 Collocation Method for Elliptic PDEs and SPDEs 

We begin by setting up Gaussian processes via reproducing kernels with differential and 
boundary operators. 

Suppose that the reproducing- kernel Hilbert space Hft-('D) is embedded into the Sobolev 
space T-^iV) where m > d/2. Let n := \m — d/2] — 1. By the Sobolev embedding the- 
orem T~L m (T>) C C n {T>). When the differential operator P and the boundary operator 
B have the orders O(P) < m — d/2 and O(B) < m — d/2, then the stochastic pro- 
cesses PS x {uj) := (Pu)(x) and BS x (uj) := (Bco)(x) are well-defined on (Qk, Fk, I 5 ' 1 )- 
According to Lemma B.l, we have PS = J2T=i(kV^kPek and BS = J2T=i CkV^kBek- 
If y, € B K (V) C rl rn (V), then P/i £ C(V) and By, £ C(dV). Lemma B.2 implies 

that P 1 P 2 K(x,y) = ET=i X l Pe k(x)Pe k (y) and B 1 B 2 K(x,y) = EZi X l Be k(x)Be k (y) 
(here we can use the fact that Cov(PS x , PS y ) = P x P y Cov{S x , S y ) and Cov(BS x , BS y ) = 
B x B y Cov(S x , S y )). Applying Lemma 2.1, we can obtain the main theorem for the construc- 
tion of Gaussian processes via reproducing kernels coupled with differential or boundary 
operators. 
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Theorem 3.1. Suppose that the reproducing kernel Hilbert space H^("D) is embedded into 
the Sobolev space T-L m (T>) with m > d/2. Further assume that the differential operator P 
and the boundary operator B have the orders 0{P) < m — d/2 and O(B) < m — d/2. 
Given a function u G Hk(T>) there exists a probability measure R^ defined on {Qk^k) = 
(~H.k(T>), £$(H.k(T)))) (as in Lemma 2.1) such that the stochastic processes PS, BS given 
by 

PS x (ta) = PS(x,u) := (Pu)(x), xeVcR d , to G Qr = Hjf(2?), 
BS w (u) = BS(x,u) := (Bu)(x), x G &D, u G U K = R K (D), 

* * 

are jointly Gaussian processes with means P\x, Bfi and covariance kernels P\P 2 K, B1B2K 
defined on (£Ik, J~k, IP^) , respectively. In particular, they can be expanded as 
00 00 
PS* = Y J (ky%Pe k {x), x G V, and BS X = J2CkVhBe k (x), x G dV, W-a.s., 

k=l k=l 

where {X k } k L 1 and {ek}^ = i are the eigenvalues and eigenf unctions of the reproducing kernel 
K and their related Fourier coefficients are the independent normal variables ( k ~ A/"(/ife, \ k ) 
and p, k := (fx, ^f\ke k ) K ,v, k G N. 

Corollary 3.2. Suppose all notations and conditions are as in Theorem 3.1. Given col- 
location points X-p := {xj}j =1 C T> and Xqx> := {a;jv+j}^i C &D, the random vector 
Spb ■= (PS X1 ,--- ,PS XN ,BS XN+1 ,--- ,BS XN+M ) T defined on (Uk^k,^) has a multi- 

normal distribution with mean rrip B and covariance matrix Kp# ; i.e., 

S PB ~ J\f(m PB ,K PB ), 
where m pB := (P//(sci), • • • , Pfi{x N ), B/j,(x N+1 ), ■■ ■ ,Bfi(x N+M )) T and 



K 



PB 



{P 1 P 2 K(x j ,x k ))^ k N =1 , (P 1 B 2 K(x j ,x N+k ))fj^ 1 
(BxP 2 k(x N+j ,x k ))^ v {B 1 B 2 K{x N+:j ,x N+k ))f^ h 



Remark 3.1. While the covariance matrix Kpp may be singular, it is always positive semi- 
definite and therefore always has a pseudo- inverse Kp B ■ 

Using Corollary 3.2, we can compute the joint probability density function p^ of Spb 
defined on (Qr, J-r, P^)- In the same way, we can also get the joint density function 
Pj of (S x ,Spb) defined on (Qr, ?Ki R^\ By Bayes' rule, we can obtain the conditional 
probability density function of the random variable S x given Spb- 

Corollary 3.3. We follow the notations of Corollary 3.2. For any fixed x G V, the random 
variable S x given Spb defined on (Hr, J~Ki P^) has a conditional probability density function 

PSW«) == %4 = -^-p (- {V - m *[ V 2 ))2 ) , v G R, • G R N+M , 
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where 

* 

m£0) := n(x) + kp B (x) T Kp B j (v - m PB ), 

a(x) 2 := K(x,x) — kpB(x) T KpB^kpB(x), 

k PB {x) := (P 2 K(x, oji), • • • , P 2 K(x, x N ), B 2 K(x, x N+1 ), ■■■ , B 2 K(x, x N+M )) . 

In particular, given the real observation y := (j/i, ■ ■ ■ ,i/n+m) T , S x conditioned on 
Spb = V has the probability density px{-\y) ■ 

This corollary is similar to the features of Gaussian conditional distributions (see [13, 
Theorem 9.9]). 

3.1 Elliptic deterministic PDEs 

Suppose that u 6 Hk-(P) is the unique solution of the deterministic elliptic PDE 

{Pu = f, inVcR d , . . 

< (3.1) 
yBu = g, on dV, 

where / : T> — >■ M and g : dT> — > R. Denote by and {yjv+fclfc^i the values of / and 

g at the collocation points Xz> and Xgx>, respectively: 

yj := f(xj), y N +k ■= g(x N+k ), j = l,-- ,N, k = 1, • • • , M. 

From now on we assume that the covariance matrix Kpp defined in Corollary 3.2 is 
nonsingular and we therefore can replace pseudo-inverses with inverses. 

Let y := (yi, ■ ■ ■ , yjy, yjv+i, ■ ■ ■ , 2//v+m) T \ and denote by Px(-\-) the conditional density 
function defined in Corollary 3.3. We approximate the solution u of (3.1) by the optimal 
estimator u(x) derived in the previous section, i.e., we maximize the conditional probability 
given the data values y : 

u(cc)»n(a;)=argmax sup Px\v\yo), x^V. 

By direct evaluation as in Section 2.2 one finds that 

u{x) := fcps(x) T Kp J B~ 1 y , x e V, 

where the basis functions kpp(x) are defined in Corollary 3.3. Moreover, the estimator 
u € Hft;('D) fits all the data values: Pu{x{) = yi, . . . , Pu(xn) = yjy and Bu(xn + \) = 
y/v+i, • • • , Bu(xn+m) = yN+M- This means that we have computed a collocation solution 
of the PDE (3.1). Also note that u can be written as a linear combination of the kernels as 
in (1.4), i.e., 

TV M 

u(x) = ^2c k P 2 K(x,x k ) +^2c N+k B 2 K(x,x N+k ), x £ V, (3.2) 
k=i k=i 
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where c := (c 1; • • • ,c N+M ) T = ^PB^Vo € R N+M . 

Finally, we can perform a weak error analysis for \u{x) — u(x)\ as in Proposition 2.2, 
and deduce that 

^ {&) = P^ (15, - u(x)\ > e\S PB = y ) = erfc {^TT 

where o~(x) 2 is defined in Corollary 3.3, and 

E e x := {u G 0, K : \u(x) - u(x)\ > e s.t. Puj(x 1 ) =y u ..., Bu>(x N+M ) = Vn+m} ■ 

Because the form of the expression for the variance c(x) 2 is analogous to that of the 
power function, we can use the same techniques as in the proofs from [11, 21] to obtain a 
formula for the order of a{x). 

Lemma 3.4. 

a(x) = 0(h™- p ~ d/2 ), xeV, 
where p := max{0(i- > ), O(B)} and hx is the fill distance of Xp and Xgx>- 

Proof Since there is at least one collocation point Xj G Xx>L)Xqx> such that ||a; — flJjHa < hx 

we can use the multivariate Taylor expansion of K(x,Xj) to introduce the order of o~(x), 
i.e., 

K(x,Xj)= -^D^D^Kix^x^x-x^ + R(x, Xj ), a,/3G< 

\a\,\/3\<n °" P ' 

where R(x, Xj) := J2\ a \,\/3\= n ■^piDfD^K(zi,z 2 )(x ~ x jT +l3 for some zi,z 2 G T> and n := 
\m — d/2] —1. The rest of the proof proceeds as in [11, Chapter 14.5] and [21, Chapters 11.3, 
16.3]. □ 

Using Lemma 3.4 we can deduce the following proposition. 

Proposition 3.5. For any e > 0, 

' im—p—d/2 x 
it i 



sup F(^)=OM , xeV, 



which indicates that 



sup P M (||u - uIIlooCu) > e ) < sup P M (££) -> 0, w/ien /ix -> 0. 

Therefore we say that the estimator u converges to the exact solution u of the PDE (3.1) 
in all probabilities P^ when hx goes to 0. 

Sometimes we know only that the solution u G T-L m (T>). In this long as the 

reproducing kernel Hilbert space is dense in the Sobolev space T-L m (T>) with respect to its 
Sobolev norm, we can still say that u converges to u in probability. 
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3.2 Elliptic stochastic PDEs 

Let £ : T> x $7^ — > R be Gaussian with mean and covariance kernel $ : P x D -> I on 
the probability space (Qyy , Tw ,^w) ■ We consider an elliptic PDE driven by a Gaussian 
additive noise £ 

fp« = / + f, inDcR rf , . . 

< (3.3) 
I Bu = g, on 9X>, 

and suppose its solution u G L^Oiyj Hft;(P)). 

Since £ is a Gaussian process, on some underlying probability space {Vtw^Fw^w) 
with a known correlation structure, we can simulate the values of £ at Xj, j = 1, . . . , N. 
Consequently, we assume that the values {yj}jLi and {yN+k\k=i defined by 

yj := f(xj) + £ Xj , yN+k ■= g(x N+k ), j = l,--- ,N, k = l,--- ,M, 

are known. 

In order to apply the general interpolation framework developed in Section 2.2, we 
consider the product space 

n KW -.= n K x n w , f kw ■- t k ® f w , f w '■= Ffl ® F w- 

We assume that the random variables defined on the original probability spaces are extended 
to random variables on the new probability space in the natural way: if random variables V\ : 
Qk P and V-i : Q\y P are defined on (fix, Fk, P^) and (Clw, Fw, Pw)> respectively, 
then 

V\(u),uj) := V\(oS), V2(oj,u>) := V^(d;), for each w G ri^- and w G fiw. 

Note that in this case the random variables have the same probability distributional prop- 
erties, and they are independent on {VI kw, Fkw-, Pyj/)- This implies that the stochastic 
processes S, PS, BS and £ can be extended to the product space {VIkw, Fkw, P^) while 
preserving the original probability distributional properties. In this case, {yi,--- ,Vn) ~ 
A/"(/,l|/), where / := (/(asi), - - - J(x N )) T and V := (*(x i9 a*))^ with * being the 
covariance kernel of £. According to Corollary 3.3 we can find the conditional probability 
density function Px{-\-), and the optimal estimator u{x) 

u{x) = argmax sup Px{ v \y s)Pw {y g) , a: G P, 

t'GK /iGH x (D) 

where 

y € :=(W,-,^) T , m 5 :=f^GR^, I := ft ^ e r^x^ 
Pw(v) : 



v / det(Z)(27r)(^+ A/ )/ 2 



expr~(«-m € ) T Et(«-m € )y v G R iV+M . 
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If Kpb + 51 is nonsingular, then one optimal solution has the form 

u{x):=kp B {x) T {K PB + Y.)- l y i , xeV, 

* 

where Kpp and kpp(x) are denned in Corollary 3.2 and 3.3. 

Similar to the analysis of the error bounds from Section 3.1, we also deduce the following 
proposition (for more details see [23]). 

Proposition 3.6. Assume ^ £ C l ^(V x V). Then 

lim sup (\\u - u\\ Lca (v) > e) = 0, for any e > 0. 

4 Numerical Experiments 

We consider the following stochastic heat equation with zero boundary condition 

' dU t = -£jU t dt + adW tti , in V := (0, 1) C R, < t < T := 1, 

C/t = 0, on <9X>, (4.1) 

J7o(a;) = i*o(x) := \/2 (sin(7rx) + sin(27rx) + sin(37rx)) , 

driven by two types of space-time white noise (colored in space) W of the form 

Wtf :=^2w t k q\(j> k , q k := — , fc (x) := >/2 sin(/c7rx), 
fc=l 

where W^, £ N, is a sequence of independent one-dimensional Brownian motions, and 
i = 1,2. Note that choosing the larger value of i corresponds to a noise that is smoother in 
space. 

The spatial covariance function R l (x,y) = X^fcLi Qk' ( i ) k(x)4>k{y)i i = 1, 2, takes the 
specific forms 

B}(x, y) = min{x, y} — xy, < x, y < 1, 

and 

i? 2 (X; y) = | + + ^ " + I^' < x < y < 1, 

l~|y 3 + s^ 3 + l x3 y - ^ x2 y + \ x y> o < y < x < i. 

The solution of SPDE (4.1) is given by (for more details see, for instance, [5]) 

oo 

Ut{x) = Y J $<t>k{x), X€V:=(0,1), 0<t<T:=l, 
fc=l 

where 



:= [ «o(x)to(s)dx, := ^ Vt + ^ / e^^d^ 

<7fc Jo 
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From this explicit solution we can get that 

r 2 

2k 2 K 2 q 2i 



E(U t (x)) = Y,&~ k2n2t Mx), Var(0i(s)) = E ^TTlgO- ~ e" 2M ) \4>k{x)\ 2 



We discretize the time interval [0, T] with n equal time steps so that 5t := T/n. We also 
choose the reproducing kernel K(x,y) := <?3,26»( x ~~ y)i where g^ l2 e is the Matern function 
with degree m := 3 and shape parameter 6 > (see Example A.l). As collocation points 
we select uniform grid points Xx> C (0, 1) and Xqx> '■= {0, 1}. Let P := I — 5td 2 /dx 2 
and B := /|{o,i}- Using our kernel-based collocation method we can perform the following 
computations to numerically estimate the sample paths ~ Ut n (xj). Algorithm to solve 
SPDE (4.1): 

1. Initialize 

• u° := (n (xi),--- ,u (x N )) T 

* / (P 1 P 2 K(x j ,x k ))fjZ v (P&Kix^XN+k))?^ 

• r\pB '■= I * * 

y(JB 1 P 2 .K'(a;jv+j ) a;fc))^= 1 , {B 1 B 2 K{x N+j ,x N+k ))fj^ li 

• B := ((P 2 K(x j ,x k ))fj!l 1 , {B 2 K{x j ,x N+k ))f k ^ 1 

• :=a 2 5t{R{x h x k ))lf =1 , Y. 

• A:= B(K PB + Z)- 1 

2. Repeat for j = 1, 2, . . . , n, i.e., for ti,t 2 , . . . , t n 

• Simulate £ ~ M (0, v|/) 

f U^ 1 + £ 




B(K PB + Z) 








Note that in the very last step the matrix A is pre-computed and can be used for all 
time steps, and for different sample paths; that makes the proposed algorithm to be quite 
efficient. 

We approximate the mean and variance of Ut(x) by sample mean and sample variance 
from s := 10000 simulated sample paths using the above algorithm, i.e., 

S S / 8 \ 

k=l i=l V k=l / 

Figure 4.1 shows that the histograms at different values of t and x resemble the theo- 
retical normal distributions. We notice a small shift in the probability distribution function 
of the solution U, at times closer to zero, and when the noise is equal to W\ (Figure 4.1, 
left panel). This shift is due to the fact that W\ is rougher in space than W 2 . 
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t = 0.325, x = 0.5082 t = 0.325, x = 0.5082 t = 0.325, x = 0.5082 t = 0.325, x = 0.5082 




Empirical - 6 - Theoretical Empirical - 6 - Theoretical 

(a) With spatial covariance R 1 (b) With spatial covariance R 2 



Figure 4.1: Empirical and theoretical probability distribution of Ut(x) for uniform points 
N := 58 and M := 2, equal time steps n := 800, 9 := 26.5, a := 1. 

Our use of an implicit time stepping scheme reduces the frequency of the white noise, 
i.e., limjf-^o SW/5t ~ So. Consequently, Figure 4.2 shows that the approximate mean is well- 
behaved but the approximate variance is a little smaller than the exact variance. According 
to Figure 4.3 we find that this numerical method is convergent as both 5t and hx are refined. 
Finally, we want to mention that the distribution of collocation points, the shape parameter, 
and the kernel itself were chosen empirically and based on the authors' experience. As 
mentioned before, more precise methods are currently not available. A rigorous investigation 
of these questions, as well as determination of precise rates of convergence is reserved for 
future work. 

5 Final Remarks 

This new numerical approach can also be used to approximate systems of elliptic PDEs 
with vector Gaussian noises ^ 1 and £ 2 or nonlinear PDEs with Gaussian noise £, i.e., 

Pu = f + £ 1 , mVcR d , fF(Pu) =V(/,0. in^cR d , 

or < 

Bu = g + £ 2 , ondV, [G{Bu) = g, on&D, 

where P := (P 1 , • • • , P n p) T is a vector differential operator and B := (B 1 , • • • , B nb ) T is a 
vector boundary operator, and F : R np R and G : R Ub — > R (see [23]). 

In addition to the additive noise case discussed here, we can also use the kernel-based 
collocation method to approximate other well-posed stochastic parabolic equations with 
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(a) With spatial covariance R 1 (b) With spatial covariance R 2 

Figure 4.2: Approximate and theoretical mean and standard deviation for uniform points 
N := 58 and M := 2, equal time steps n := 800, 6 : = 26.5, a := 1. 



multiplicative noise, e.g., 

dU t = AU t dt + ip(U t )dW t , in£>cM d , < t < T, 
BU t = 0, on &D, (5.1) 

^U = u , 

where : R -»• R. Since / t V_ ^(E7,)dW, « ^(Ut^SWj, the algorithm for SPDE (5.1) is 
similar to before: 



1. Initialize 



u° := (u (xi), • • • ,u {xn)) 



T 



P B ' — I 5^ 

j5iP 2 -^(a;Af+i 5 a; fc))^fcfi 5 { B i B 2K{x N+j ,x N+k ))f^ l/ 
• B:= ((^xj))^, (Ba^(*i,^+*))JJiS) 

2. Repeat for j = 1, 2, • • • , n, i.e., for t±, i 2 , ■ ■ ■ , t n = T 

. V! := diag faut 1 ), ■ ■ ■ M^N 1 )) , V 2 := °j 
. V|/ : = V1V0V1, E := V 2 Z V 2 
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Figure 4.3: Convergence of mean and variance with respect to refinement of points and 
time steps for a := 1. (The relative RMSE of exact U and approximate U is defined by 



RMSE(C7, U) := £Li " Ufa, x k ))*/\Mt kl -)ll§o-) 



• Simulate £ ~ N (0, V) , A := B(K PB + Z)" 1 

Of course, now the matrix A needs to be updated for each time step and for each sample 
path so that the algorithm is much costlier. 

A Reproducing-Kernel Hilbert Spaces 

Definition A.l ([21, Definition 10.1]). A Hilbert space B. K (V) of functions / : V ->■ R is 
called a reproducing-kernel Hilbert space with a reproducing kernel K : 2? x D — > M. if 

(i) lT(.,y) G Hjf (X>) and (u) /(y) = (if (•, y), /) A ^, 

for all / € H^'('D) and all y 6 "D. Here (•, -)x,x> is the inner product of H^('D). 

According to [21, Theorem 10.4] all reproducing kernels are positive semi-definite. [21, 
Theorem 10.10] shows that a symmetric positive semi-definite kernel K is always a repro- 
ducing kernel of a reproducing-kernel Hilbert space H;c(D). 

If D is open and bounded (pre-compact) and K £ L2(T> x T>) is symmetric positive 
definite, then Mercer's theorem [11, Theorem 13.5] guarantees the existence of a countable 
set of positive values Ai > A2 > • • • > with YlT=i Afc < 00 and an orthonormal base {efc}^ =1 
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of Ij2(T>) such that K possesses the absolutely and uniformly convergent representation 

oo 

K(x, y) = Xkek(x)e k (y), x, y G V. 
k=l 

This Mercer series of K implies that 

Afce fc (t/) = / K(x,y)e k (x)dx, y G V, k 6 N. 
Jv 

Definition A. 2. The sequences {X k } k x L 1 and {e k } k x L 1 given above are called the eigenvalues 
and eigenf unctions of the reproducing kernel K. 

Since {e k } k x L 1 is orthonormal in L2(D) we can compute the series expansion of the 
integral-type kernel K defined in Lemma 2.1, i.e., 

oo oo „ oo 

K(x,y) = ^2Y^ / ^je j {x)e j {z)X k e k (z)e k (y)dz = ^ \ 2 k e k {x)e k {y). 
j=i k =i J v k =i 

It is easy to check that K S L2(£> x T>) is symmetric positive definite and { A^}^ =1 and 
{e k } k x L 1 are the eigenvalues and eigenfunctions of K. 

Theorem A.l ([21, Theorem 10.29]). Suppose that K 6 xT>) is a symmetric positive 
definite kernel on a pre- compact T> C K d . Then its reproducing-kernel Hilbert space is given 
by 





2 1 


/ f(x)e k (x)dx 


< oo 


Jv 





and the inner product has the representation 

(f,g)K,v = y^T- / f(x)e k {x)dx I g(x)e k (x)dx, f,g£ Hff(P), 
fc=i fc 

where {Afc}^l 1 and {e/ c }^ =1 are i/ie eigenvalues and eigenfunctions of K. 

We can verify that {^/Xkek}^^ is an orthonormal base of Hk(T>). 

Example A.l. The papers [9, 22] show that the Sobolev spline (Matern function) of degree 
m > \ , 

ol— m— d/2 

9m,eix) ■■= ^r^)^-^ ^^^)^ 72 ^/^^^^). *6R< 8>0, 

is a full-space Green function of the differential operator L := (6 2 I — A) m , i.e., Lg m fi = 5q, 
where 1 1— > K y (t) is the modified Bessel function of the second kind of order v. The kernel 
function 

K m 6 (x,y) := g m ,e(x - y), x,y eR d , 
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is a positive definite kernel whose reproducing-kernel Hilbert space is equivalent to the L2- 
based Sobolev space of degree m, i.e., H^ mfl (R d ) = Ti m (R d ). Its inner product has the 
explicit form 



</,<?} 



[ P e f(x) T P g g(x)dx, f,geR Kmg (R d ), 



where := (Q°, Q 1 , , Q m ) and 



Q J -={ A , T , a i : =\h? 7T' k ^ n o, J = 0,l,---,m. 

{ aj A k V T , j = 2k + l, \3 ] -{m-jy- 

According to [3, Theorem 1.4.6], the reproducing-kernel Hilbert space Hjr mei (X>) is endowed 
with the reproducing-kernel norm 

\\f\\ Km , e ,v = _ mf {\\f\\ Kme ^ ■ f\v = A , / € H Am9 (P). 

Moreover, if the open bounded domain 2? C M. d is regular then Hjf mfl (I?) is equivalent to 
the L2-based Sobolev space of degree m, i.e., H^ me (P) = % m (T>). 



B Differential and Boundary Operators 

In this section we define differential and boundary operators on Sobolev spaces T-L m (T>). The 
differential and boundary operators in this paper are well-defined since we assume that the 
open and bounded domain D is regular, i.e., it satisfies a strong local Lipschitz condition 
which implies a uniform cone condition (see [1, Chapter 4.1]). This means that T> has a 
regular boundary &D. 

Let the notation for typical derivatives be 

k=l k k=l 

We extend these derivatives to weak derivatives (see [1, Chapter 1.5]) using the same symbol 
D a . Using these weak derivatives, the classical L2-based Sobolev space T-L m (V) is given by 

H m (V) : = {/ G L[ oc (V) : D a f G L 2 (V), \a\ < m, a G N^} , m G N , 

equipped with the natural inner product 

(f,g) m ,v:= V / D a f(x)D a g(x)dx, f,g£H m (V). 

|a|<m 

The weak derivative D a is a bounded linear operator from T-L m {T>) into I^I*) when |a| < m. 
Moreover, D a f is well-posed on the boundary dT> when / G C m (T>) and \a\ < m — 1 and we 
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denote that D a \ox>f := D a f\gx>- The book [1] and the paper [10] show that D a \ox> can be 
extended to a bounded linear operator from 7-L m (T)) into L2(9X?) when \a\ < m — 1 because 
T> has a regular boundary dT>. The L2 (<9X?)-inner product is denoted by 

{fid)&D= / f( x )g( x )dS(x), when d > 2 and dT> is the boundary manifold, 
Jdv 

and 

(/, g)dv = f(b)g{b) + f(a)g(a), when d = 1 and dV = {a, b}. 
Definition B.l. A differential operator P : H m (T>) — > L2(T>) is well-defined by 

P= c ^ Q > where c « G C°°(V) and q£< m G N 0j 

\a\<m 

and its order is given by O(P) := max ||a| : c a ^ 0, \a\ < m, a G Nq} . A boundary oper- 
ator B : H m (V) -> L 2 (d£>) is well-defined by 

5 = ^ b a D a \ dv , where 6 Q G C°°(d£>) and a G Nq, m G N, 

|a|<m— 1 

and its order is given by O(B) := max { |a| : 6 a ^ 0, |a| < m — 1, a G Nq} . 

It is obvious that the differential operator P and the boundary operator B are bounded 
(continuous) linear operators on T-L m {T>) with values in L2 whenever 0{P) < m and O(B) < 
m — 1. Much more detail on differential and boundary operators can be found in [1, 10]. 

If Hr-(P) is embedded 1 into T-L m {V) then the eigenvalues {A/ C }^1 1 and eigenfunctions 
i e k}kLi °f the reproducing kernel K satisfy 

^fcll e fcllm,z> - C 2 \\v^k~ e k\\K,v = C 2 , k G N, 

because {\/Afc e fc}fcLi is an orthonormal base of H^(X?). When m > d/2 then H m (T)) is 
embedded into 0(2?) by the Sobolev embedding theorem. This implies that K G C(T>xD) c 
L^l? x X?) because K(-,y) G C(X?) for each y £ T> and is symmetric. Based on these 
properties, we can introduce the following lemma. 

Lemma B.l. Consider a differential operator P with order 0(P) < m and a boundary 
operator B with order O(B) < m — 1, where m > d/2. If the reproducing-kernel Hilbert 
space Hftf('D) is embedded into the Sobolev space H m (V), then 

00 00 
Pf = Y,f^kPe k , Bf = J2fkV%Be k , f G H K (V), 

k=t k=l 

where fk = {f-,&k)K,v for each k G N and {A/ c }^ =1 and {e^^i are the eigenvalues and 
eigenfunctions of the reproducing kernel K. 

The recent papers [9, 10, 22] show that there exist kernels K whose reproducing-kernel Hilbert space 
Hr-(2?) is continuously embedded into the Sobolev space H m (X>). 
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Proof. According to Theorem A.l each / £ Hr-(£>) can be expanded as / = YlkLi fk\/Xk e k- 

* 2 

Since {\/Afcefc}^ =1 is an orthonormal basis we have /fc < °°- Let /„ := 53fc=i fkV^~k e k 

for each n 6 N. Then 



2 



0, when n — > oo. 



||/n-/||^<^ 2 H/n-/ll^<C 2 £ |/ fc 

k=n+l 

The proof is completed by remembering that P and B are bounded linear operators on 
H m {V). □ 

If Hjf(P) is embedded into H m {T>), then for each \a\ < m, < m and a,/3e Nq, we 
have 

1/2 



fe=i 



da;dy 



<E ^ll £>Q e fc || L2(2?) ||^e fc || L2(23) < £Ai||e fc ||^ < C 2 f> < oo, 

k=l k=l k=l 

which implies that K G n m > m {V x £>). Let n := [m - d/2] - 1. The Sobolev embedding 
theorem shows that H m > m (T) xD) C C n,n (P x P). Then we can obtain the following lemma. 

Lemma B.2. Consider a differential operator P with order O(P) < m—d/2 and a boundary 
operator B with order O(B) < m — d/2, where m > d/2. If the reproducing-kernel Hilbert 
space H/^(D) is embedded into the Sobolev space % m (P), then 



oo 



P 1 P 2 K(x,y) := P Zl P Z2 K( Zl ,z 2 )\ Zl=x , Z2=y = ^XlPe k (x)Pe k (y), 

k=i 

oo 

* * 

BiB 2 K(x,y) := B Zx B Z2 K{z 1 , z 2 )\ Zl=x>Z2 = y = ^ X 2 k Be k (x)Be k (y), 

k=i 

oo 

P 1 B 2 K{x,y) := P Zl B Z2 K(z!, z 2 )\ Zl=x , Z2=y = ^ X 2 k Pe k (x)Be k (y), 



k=i 



where {X k } k x L 1 and {e k } k x L 1 are the eigenvalues and eigenf unctions of K . Moreover, P\P 2 K £ 
C(P x V), B\B 2 K G C(dV x dV) and P X B 2 K G C(P x SP). 
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